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^\ ' The stability and melting transition of a single layer and a bilayer crystal consisting of charged 

0^ . particles interacting through a Coulomb or a screened Coulomb potential is studied using the Monte- 

0^ ' Carlo technique. A new melting criterion is formulated which we show to be universal for bilayer 

as well as for single layer crystals in the case of (screened) Coulomb, Lennard-Jones and 
repulsive inter-particle interactions. The melting temperature for the five different lattice structures 
^ ' of the bilayer Wigner crystal is obtained, and a phase diagram is constructed as a function of the 

I interlayer distance. We found the surprising result that the square lattice has a substantial larger 

melting temperature as compared to the other lattice structures. This is a consequence of the 
specific topology of the defects which are created with increasing temperature and which have a 
larger energy as compared to the defects in e.g. a hexagonal lattice. 
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I. INTRODUCTION 

Phase transitions and in particular the melting transition is one of the most fundamental problems in condensed 
matter and statistical physics. The most intensively studied model system is the one consisting of charged particles 
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(— ( ; with a pure spherical symmetric repulsive inter-particle interaction (see e.g. Examples are the electron Wigner 
Q ■ solid, colloidal spheres and dusty plasmas. After the discovery of the Wigner crystallisation of electrons above the 
^ O ^ , surface of liquid helium by Grimes and Adams Q the interest in the melting of low dimensional systems has increased. 

Recently the layered close-packed crystalline structure and their structural transition were observed in experiments 
T— I , with dust rf discharge [|j and with ion layered crystals in ion traps Q . Motivated by the theoretical works of Nelson, 
K*" ' Halperin ||], and Young |^ who developed a theory for a two stage continuous melting of a two dimensional (2D) 
crystal and which was based on ideas of Berenzinskii [Q , Kosterlitz and Thouless , several experimental p|~p^ and 
' theoretical studies p"^"|25[| were devoted to the melting transition of mainly single layer crystals. In this case the 
QQ \ hexagonal lattice is the most energetically favored structure for potentials of the form 1 / r" [ p6p^ ] . The effect of the 
' dimensionality of the system on the melting transition is another important issue which has been investigated recently. 
On An important subsystem in this field is the bilayer system consisting of two parallel layers of charged particles with 
a pure repulsive inter-particle interaction. The interlayer distance is an additional variable with which the effective 
interparticle interaction can be altered. 

While the classical single layer can crystallise only in a hexagonal lattice it was found that the bilayer system 
exhibits a much richer structural phase diagram. Previously, the different types of lattices and structural transitions 
in a multilayer crystal at zero temperature with parabolic lateral confinement was analysed in [p8|-^ . Different classes 
^ of lattices of the double-layer crystal were specified in |3l) as function of the interlayer spacing and electron density. 
O ■ A systematic study of the stability and the phonon spectra of these crystal phases was presented in Ref. . The 
, collective modes of the corresponding bilayer liquid state was recently studied in p3] . The extension of the classical 



bilayer model to the finite vertical coupled classical dot system was presented in Ref. |3J] where a rich collection of 
structural first and second order phase transitions was found. The quantum bilayer system was investigated by several 
5^ groups Hip. 

5^ , Here we will concentrate on the melting of the classical bilayer crystal and study the influence of the different crystal 
structures on the melting temperature. The equilibrium states of a bilayer crystal at different temperatures (T) will 
be investigated using the Monte Carlo (MC) simulation technique. Preliminary resuls of this study were published in 

ii- . . . . . ^ . . 

The melting of the classical bilayer crystal was studied in Ref. [ p2[ using the Lindemann criterion where the 
mean square displacement was calculated within the harmonic approximation. It is well-known that the harmonic 
approximation can only give a rough estimate of the melting temperature. For example, for the single electron layer 
using the Kosterlitz-Thouless-Halperin-Nelson- Young (KTHNY) theory and inserting the Lame coefficients of the 
T = lattice one obtains F = 78, while inclusion of non linear effects, using Monte-Carlo simulations |l^], results into 
F = 125 which compares to the experimental result F = 137 ± 15 (see e.g. [Eoi). This is the motivation to go beyond 
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the harmonic approximation and to investigate the importance of non hnear effects on the melting temperature. Using 
the Monte Carlo simulation technique has the added advantage that information can be obtained about the topology 
of the defects which are responsible for melting. 

The present paper is organised as follows. The model system and the numerical approach is given in Sec. II. In Sec. 
Ill we construct the solid-liquid phase diagram for Coulomb bilayer crystals. In Sec. IV the influence of screening 
on melting is studied and we formulate a new criterion for the melting temperature. In Sec. V we show that it gives 
accurate results for the melting temperature when we use the new melting criterion and the harmonic approximation. 
The existence of hysteresis during melting of the 2D single layer crystal is considered in Sec. VI. The topology of 
defects in the hexagonal type crystal is considered in Sec. VII. The gallery of defects for the square bilayer crystal is 
presented in Sec. VIII. Our results are summarised in Sec. IX. 

II. THE MODEL SYSTEM AND NUMERICAL APPROACH 

In the present study we limit ourselves to infinitely thin bilayers with equal density, n/2, of charged particles in 
the two layers. In the crystal phase the particles are arranged in two parallel layers in the (x, ?/)~plane which are a 
distance d apart in the z-direction. The lattice structure in the top and bottom layers are the same, but the particle 
sites in opposite layers are staggered (i.e. closed-packing system). A single layer crystal is a limiting case of a bilayer 
crystal with d = and particle density n. 

We assume that the particles interact through an isotropic Coulomb or screened repulsive potential 

" , I ^ I exp(-K| fi - fj I), (1) 

where q is the particle charge, e is the dielectric constant of the medium the particles are moving in, r — {x,y, z) the 
position of the particle with r —\ r |, and 1/k is the screening length. The type of lattice symmetry at T = depends 
on the dimensionless parameter v = d/ao, where oq — l/^/im is the mean interparticle distance. For the classical 
Coulomb system (k = 0) at T ^ there are two dimensionless parameters u and T = q^ / ea^kBT which determine the 
state of the system. The classical Yukawa system (k > 0) at T ^ is characterised by three independent dimensionless 
parameters: v, T and A = kcq. Below we measure the temperature in units of Tq = (f' /^o-okB and the energy in 
Eq = ksTo. In addition, we will also consider the melting of a single layer crystal with a Lennard-Jones (l/r^^ — 
inter-particle potential and with a 1 /r^^ repulsive potential as well. 

In our simulations the density of particles n remains the same for all interlayer distances v while the interlayer 
distance changes. The initial symmetry of the structure is set by the primitive vectors, the values of which are derived 
from a calculation of the minimal energy configuration for fixed at T = using the Ewald method |4^ ]. In our 
numerical calculations we took a fragment of lattice (ranging from N ~ 228 to = 780 particles) and used periodic 
boundary conditions. 

The structure and potential energy of the system at T are found by the standard Metropolis algorithm in 
which at some temperature the next simulation state of the system is obtained by a random displacement of one of the 
particles. If the new configuration has a smaller energy it is accepted and if the new energy is larger the configuration 
is accepted if <5 > r, where the probability S = exp{— AE / ksT) with AE the increment in the energy and r is a 
random number between and 1. We allow the system to approach its equilibrium state at some temperature T, 
after executing (10'* 5 x 10^) 'MC steps'. The potential energy at T 7^ is found by summation over all particles 
and their periodic images using the Ewald method . 

III. SOLID-LIQUID PHASE DIAGRAM 

In it was found that the bilayer Coulomb crystal at T = exhibits five different lattices as function of the 
interlayer distance: < 0.006-hexagonal, 0.006 < ly < 0.262-rectangular, 0.262 < ly < 0.621-square, 0.621 < ly < 
0.732-rhombic, and 1/ > 0.732-hexagonal. 

It should be noted that the behaviour of the system near the melting point is essentially non-linear and therefore 
the melting temperature obtained in Ref. |Q has only qualitative predictive power. Therefore, we revise the melting 
phase diagram and present a more accurate calculation of it, by performing MC simulation of the melting transition. 

We use different criteria to find the critical melting temperature. The parameter =< u1 > /a§, where < ul > is 
the mean square displacement of the particles from their equilibrium positions was introduced in Esl to characterise 
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the system order. It is well-known that L* diverges for a 2D system. We will use a modified parameter L as introduced 
in Ref. |p^ , which enables us to minimize the system size efi^ect which gives significant contributions in the usual 
parameter L*. For the modified parameter L —< v? > /aq, < > is defined by the difi^erence in the mean square 
displacement of neighboring particles from their initial sites tq and can be written as 

< (nJ2j^J2((^^^ - - (^^- - ^j-o))')' (2) 

i=l " j=l 

where <> means averaging over MC steps, the index j denotes the Nnb nearest neighbours of particle i in the upper or 
in the bottom layers. For the hexagonal lattice Nnb — 6, while for rectangular, square and rhombic lattices Nnb — 4. 

For sufficiently low temperatures the particles exibit harmonic oscillations around their T = equilibrium position, 
and the oscillation amplitude increases linearly with temperature, which leads to a linear increase of L with T. Figs. 1- 
3(a) show the modifies parameter L as a function of the reduced temperature T/Tq: i) for the hexagonal single layer 
crystal {v = 0) (see Fig. 1(a)), ii) for the bilayer crystal with a square lattice {i/ = 0.4) (see Fig. 2(a)), and iii) for 
a hexagonal type bilayer crystal (i/ = 0.8) (see Fig. 3(a)). At higher temperatures non-linear effects are important 
and L becomes a nonlinear function of T. The particle oscillation amplitude increases faster than linear with T, but 
the system is not melted, since the particles display still an ordered structure. The rapid increase of L with T is 
a manifestation of the anharmonicity of the motion of the particles. Melting occurs when L increases very sharply 
with T. In our previous study | |3^ we assumed that melting occurs when L ~ 0.1 as was obtained from numerical 
simulations of melting in a single layer crystal |l6[ . From the present study we find that L ~ 0.1 is not a good criterion 
for melting, except for the case ly — while for ly = 0.4 and — 0.8 the system is still ordered when L ^ 0.1. Note, 
that we defined L —< v? > /aq, where ao is the mean inter-particle distance for a single layer crystal with density n, 
for all interlayer distances v. In the case when j/ — s- oo and the inter-layer particle interaction is negligible small, we 
obtain L ~ 0.2 at melting because now the particle density in each layer is y/2 times smaller. 

As a second independent parameter from which we determine the melting temperature we calculated the height of 
the first peak in the pair correlation function as function of temperature. The pair correlation function is calculated 
for the particles in one of the layers (see Fig. 4). This is shown in Fig. 4 for an interlayer separation v = 0.4 and 
three different temperatures. The height of the first peak in g(r) is plotted as function of temperature in Figs. l-3(b) 
for v = 0, 0.4, 0.8, respectively. The value of the first peak of the pair correlation function decreases with increasing 
temperature, and exhibits a quick decrease at T = Tme/- This is most pronounced for the case of a single layer 
crystal (Fig. 1(b)). Even in the liquid phase the pair correlation function still has a peak structure due to the local 
inter-particle correlation. In Fig 4 the state of the system just before melting is denoted by the doted curve and the 
dashed curve refers to the state after melting. 

In the crystalline state the potential energy of the system increases almost linearly with temperature and only a 
weak non linear T-dependence is observed. At larger temperature a quick rise in the energy is clearly seen in Figs. 1- 
3(c) denoting the beginning of melting and is related to the unbinding of dislocation pairs, which we will discuss 
below. At melting the potential energy has a very steep rise. The heat capacity of a 2D structure is determined 
by C = {de/dT)v=const, where the total energy e = ksTo + E includes the kinetic and potential energy. The heat 
capacity of the crystal and the liquid phases near melting are very close to each other in two dimensions, and the 
first derivative of the energy with respect to T taken in the crystal phase and the liquid phase are practically the 
same. The potential energy decreases with increasing u which is a consequence of the decreasing inter-layer Coulomb 
interaction. Note, that the potential energy difference between the liquid and solid phases near melting for the square 
bilayer crystal dX v — 0.4 is of size (5e = 0.71 x lO^^fcsTo which is about a factor 2 larger than for a hexagonal 
lattice, i.e. v — 0, 5e = 0.39 x 10~^A:Bro, and ai u = 0.8, 5e = 0.31 x 10~^A:Bro- Moreover, as follows from the 
temperature behaviour of the parameter L, the peak of the pair correlation function and the potential energy, the 
square lattice turns out to exibit a substantial higher melting temperature, and consequently is more stable against 
thermal fluctuations, than the hexagonal one. The enhanced stability of the square lattice can be understood from 
an analysis of the defects which are formed during melting and which will be discussed in Sees. VII and VIII. 

In order to determine the type of melting transition and the melting point the bond-orientational and translational 
correlation functions are calculated ||l^,|l^,^ which sometimes can also be obtained experimentally |^,|l2|,^. But 
our finite fragment of 780 particles (with periodic boundary conditions) is too small for a reliable analysis of the 
asymptotic decay of the density correlation functions. Therefore, we calculate the bond-angular order factor Go 
which originally was introduced by Halperin and Nelson [Q. For each layer {k) we define the bond-angular order 
factor 
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^0 = (]vE]v^E exp(*iV„,0,,„)), (3) 

j=l " n=l 

and the translational order factor 

2 

where 0j,„ is the angle between some fixed axis and the vector which connects the jth particle and its nearest nth 
neighbour, and G is a reciprocal-lattice vector. The latter we took equal to the smallest reciprocal lattice vector. The 
total bond-angular order factor of the bilayer crystal is defined as Gg — (Gg + Gg)/2 and similar for Gtr- 

These order factors are shown in Figs. 1-3 (d) as function of temperature. For low temperatures the factors Gg 
(solid circles) and Gtr (open circles) decrease linearly with T, and when Gg reaches 0.45, both order factors drop 
quickly to zero at the same temperature as seen in Figs. l-3(d). From the behaviour of the order factors we can 
derive the temperature at which order is lost in the system. Our numerical results show that the bond-order angular 
order factor: 1) decreases linearly with increasing temperature, except very close to the melting temperature, where 
it decreases faster, and 2) it drops to zero just after it reaches the value of 0.45. 

In the phase diagram of Fig. 5 we show the melting temperature as a function of ly, which is the interlayer distance 
for fixed particle density. All criteria mentioned above gave the same temperature Tmei- For v = we obtained 
the critical F* = 132, resulting in T/Tq = 0.0076. This critical value was first measured in Ref. ||] and found to 
be 137 ± 15. Later experiments ^] and simulations [Q gave critical F* within this range. As seen in Fig. 5 the 
hexagonal (I and V), the rectangular (II) and rhombic (IV) lattices melt at almost the same temperature. This implies 
that the inter-layer correlation does not strongly influence the melting temperature, although it determines which 
lattice structure is stable and has the lowest energy. Note, however that the melting temperature for phase IV and V 
are slightly smaller than the corresponding temperature for phases I and II in the given ly range. We checked that in 
the hmit T^ei/To ~ 0.0076/^2 « 0.0537, which is reached for = 2 ^ 3. 

The general shape of the phase diagramm agrees with the one obtained earlier using the Lindemann criterion 
where the mean square displacement was calculated within the harmonic approximation. But there are several 
important differences: 1) the absolute value of the melting temperature was underestimated in Ref. p2[ by about a 
factor of 2; 2) at the second-order structural phase transition points the melting temperature becomes zero wich was a 
consequence of the softening of a phonon mode at T = 0. With the present Monte-Carlo simulations we find Tmei 7^ 
which is due to the importance of non-linear effects; 3) in Ref. the maximum melting temperature in phase II was 
larger than the one in phase III (square lattice) which is opposite to what is found in the present calculation; and 4) 
the melting temperature in phase IV was smaller than in phase V for v < 1, which is opposite to the present results. 

For the square bilayer crystal (III) the melting temperature increases with rising v and only for v > 0.4 we 
found that Tmei starts to decrease with increasing v. It is surprising that the square lattice bilayer crystal has the 
maximal melting temperature Tmei = 0.01078To which exceeds the critical temperature of the single layer crystal 
Tmei — 0.0076To, where the particle density is two times larger. The most stable square lattice bilayer occurs at 
v = 0.4 which results into the critical coupling parameter F* = 187. In the present paper we will not discuss the 
temperature induced structural phase transition, but limit ourselves to the melting transition. We denote in Fig. 5 
the phase boundaries between the different crystal phases by vertical dotted lines. The analysis of the stability of 
the classical bilayer crystal was carried out in Ref. |^2| using the harmonic approximation and it was found that 
near the structural phase transition the melting temperature was strongly reduced due to the softening of a phonon 
mode which may lead to a re-entrant behaviour. No such behaviour was found in the present Monte-Carlo simulaion. 
Near the structural phase transition we found deformed lattice structures with a long wavelength deformation. At 
present we cannot draw any definite conclusions from it and we need to simulate larger crystal fragments during more 
Monte-Carlo steps. But this observation docs not influence our conclusions concerning the melting behaviour. 

IV. MELTING OF BILAYER CRYSTAL WITH SCREENED INTER-PARTICLE INTERACTION AND 

FORMULATION OF A NEW MELTING CRITERION. 

Symmetry of the lattice of the bilayer crystal at T = might depend not only on the interlayer distance, but also 
on the type of the inter-particle interaction. The phase diagram corresponding to the lowest energy lattice structures 
for a bilayer crystal at T = with screened Coulomb inter-particle interaction is shown in Fig. 6 as function of the 
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screening length A. Notice that the phase boundaries depend only weakly on the screening length. On the other 
hand, the melting temperature reduces considerably with screening as is apparent from Fig. 5 in the case A = 1 (solid 
circles) and A = 3 (open triangles). The melting temperature of a single layer crystal with Coulomb (i.e. A = 0) 
interaction is T^ei = 0.0076ro, with screened one is T^ei = 0.0066ro for A = 1, and for A = 3 it is T,nei = 0.0035To. 
Note, that the bilayer crystal with screened inter-particle interaction has the same qualitative melting curve as the 
Coulomb bilayer crystal (see Fig. 5). 

For the case of a screened Coulomb interparticle interaction (A = 1) the change of the potential energy with 
temperature as referred to the T = result is shown in Fig. 7(a) for two different inter-layer distances. Notice 
that the potential energy initially increases linearly with temperature but then near the melting temperature it rises 
quickly within a very narrow temperatures range. It is apparent that the square lattice bilayer crystal (v = 0.4) 
has an enhanced melting temperature as compared to the hexagonal lattice (i.e. v = 0). The system looses order 
just after the bond-angular order factor reaches the value 0.45 as is apparent from Fig. 7(b). Also, for the screened 
Coulomb interaction the square bilayer lattice exibites an enhanced melting temperature. 

From the present numerical results for Gg for a screened Coulomb interaction and from previous section for a 
pure Coulomb interaction for single and bilayer lattices we formulate a new criterion for melting which we believe is 
universal: melting occurs when the bond-angle correlation factor becomes Gg ~ 0.45. We also checked the validity of 
this criterion for a single layer crystal with a Lennard- Jones V — l/r^^ — 1/r^ and a repulsive V = interaction 
potential. The results for the translational and the bond-angular factors for single and bilayer systems and for 
different inter-particle interactions are shown in Fig. 8 as function of the ratio T/Tmei- These results confirm the 
universality of the proposed criterion. It is valid for single and bilayer systems and we believe it is also independent 
of the functional form of the inter-particle interaction. 



V. HARMONIC APPROXIMATION 



Given the fact that the bond-order angular order factor decreases almost linearly with increasing temperature, 
and melting occurs when the bond-angular order factor equals 0.45 we will calculate the melting temperature within 
the harmonic approximation. Therefore, we reformulate the definition of the angular-bond order factor within the 
harmonic approximation in terms of eigenfrequencies and eigenvectors of the T — Q phonon spectrum. Previously 
such a method was used in Ref. in the calculation of the modified Lindcmann parameter for a finite cluster and 
in Ref. for the bilayer system. In the present paper we apply the same algorithm in the calculation of the bond 
angular order parameter. We consider a finite crystal fragment with periodic boundary conditions and diagonalize 
numerically the corresponding Hessian matrix. The theoretical background can be shortly described as follows. Any 
thermodynamically averaged observable which is a functional of the particles coordinates r; can be written as 

G^ j G{R)exp{-U{R)/T)dR, 

where R — (ri, r2, •••7 r^f) is a multidimensional vector describing the system coordinates and U is the potential energy. 
Within a standard linear approach we expand both G(i?) and U (R) near the equilibrium state R = Rq + ^, where 
dU / dR\j^^j^^ = 0, and obtain 

The multidimensional integration can be performed after diagonalisation of the Hessian matrix which has the elements 

d^G 

dm,' 

and after introduction of the normal mode coordinates. 

Thus within the harmonic approach the angular-bond order factor is given by the following expression 

Ge = l- ^ XI ((^" " ^"i) ^ ^ ~ ^ ^'=1) ^ ~ ^fei.*)) (5) 

n=l fc=l 1=4 

where tOi, Ai are the eigenfrequencies and the eigenvectors of the Hessian matrix formed by the second derivative of 
the potential energy; the indices nl, kl refer to the nearest neighbours of the particles n, fc, respectively, where the 



5 



summation is performed over all particles. The first three eigenmodes have uji = and correspond to the center-of- 
mass motion and rotation of the system as a whole which we removed from the summation in Eq. (5). The melting 
temperature is then obtained by taking Gg equal to the value 0.45. The obtained results are shown in Fig. 9 by the 
solid curves for the square lattice (Fig. 9(a)) as function of the interlayer distance and the hexagonal (Fig. 9(b)) lattice 
{u = 0) as function of the screening length. Thus using our new criterion together with the harmonic approximation 
reproduces Tme; rather well and agrees with our MC calculations within 10% for 0.35 < ly < 0.55 in the case of the 
square lattice. Near the structural transition the harmonic approximation underestimates the melting temperature 
due to the softening of the T — phonon modes (see Ref . ) . 

VI. IS THERE HYSTERESIS AT MELTING? 

One of the most interesting questions in the field of phase transitions is whether or not there is a hysteresis 
behaviour at melting and freezing of the crystal structure. Even if such an hysteresis is found it can be an artefact of 
the numerical simulation arising from the fact that the equilibrium state of the system is not reached. To shed light on 
this problem we studied more carefully the melting of the single layer crystal in the temperature range where a rapid 
increase of the potential energy was found Ti < T kT^, (see Fig. 1(c)). We took very small temperature increments 
of O.OOlTmej- At the temperature Ti — 0.007560To, which is just below the melting point, during numerical annealing 
of the system which lasts for 5 x lO^iV MC steps, the potential energy practically remains the same and changes only 
with 3% of the total energy jump 5e- Note, that the Monte-Carlo 'time' step is about a factor 10 larger than that in 
typical molecular dynamic calculations. At the melting temperature = 0.007565To the energy rises rapidly during 
10^ MC steps, and then the system remains in the liquid state during 5 x 10^ MC steps. We also found such a state of 
the system at the intermediate temperature T2 = 0.007562To when the system switches from solid to liquid and back 
as function of the simulation time. The system exhibits a network of grain boundary defects which makes it possible 
for the system to get back to a defect- free structure at the same temperature (see Fig. 10). Such a type of equilibrium 
behavior of the system at the critical point was found earlier by Morf . The potential energy of the system at 
this temperature and the order factors Ge, Gtr are shown in Fig. 11 as function of the 'MC time'. As function of 
time the system switches between the solid and liquid state where it looses and regains translational and orientational 
order. The existence of this state shows that there is a temperature range (at least for a crystalline fragment of 448 
particles) in which the crystalline and liquid states have equal probability to occur and, consequently, no hysteresis 
occurs in the Coulomb single layer crystal. It should be noted that the potential energy changes very rapidly, in an 
extreme narrow range of temperatures. Whether or not this rapid rise is continuous or discontinuous is not yet clear. 

VII. DEFECTS IN THE HEXAGONAL LATTICE. 

In this section we consider the topology of the defects which are created during melting of a single layer crystal. 
At low temperature the probability of defect creation is small and the lattice remains perfect. With rising T the 
energy increases and the system can overpass some potential barrier and we found that during the MC simulation 
the system transits from one metastable state to another. The different metastable states differ by the appearance 
of isomers in the crystal structure which appear with different probabilities. In order to find these isomers in the 
equilibrium state of our system we choose, during our MC simulation, a number of instant particle configurations. 
Then we quickly freeze each instant particle configuration which enables us to exclude the 'visual' oscillation effect, 
i.e. the displacements of the particles from their lattice configurations due to thermal fluctuations. After freezing 
each configuration we find the topology and energy of the isomer and the bond-angular and the translational order 
factors of this crystal structure. 

The isomers at Ti = 0.00756ro, which is a temperature just before melting, at — 0.00765ro which is just above 
melting, and at the intermediate temperature T2 — 0.007562 (see Fig. 1(c)), were obtained. Each point in Fig. 12(a) 
represents one configuration containing an isomer in a single layer (y — 0). Similar results are shown in Fig. 12(b) for 
the square lattice bilayer system {v = 0.4). 

With rising temperature first, point defects such as a 'centered vacancy' and a 'centered interstitial' are created 
at T < Ti — 0.00756To. These calculated defects as obtained by freezing the instant particle configuration to zero 
temperature are shown in Fig. 13(a,b). Point defects appear in pairs in our MC calculations, which are a consequence 
of the periodic boundary condition. The total energy of a non bounded pair consisting of a 'centered interstitial' and 
a 'centered vacancy' is E = 0.29kBTo. We checked the accuracy of this value on fragments of 512 and 780 particles 
which contain only one pair of 'centered interstitial' and 'centered vacancy'. In both cases the minimum energy of a 
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pair of uncorrelated point defects was the same within 1%. The quartet of bound dischnations (two sevenfold and 
two fivefold) also appear at this temperature (see Fig. 5(a) of Ref. |^) and has an energy E = 0.285A:_bTo. The point 
defects and dischnations bound into a quartet have only negligible effect on Gg and Gtr (see Fig. 12(a), group 1). 
The order factors of the 'cold' configurations (freezed to T = 0) containing a 'centered interstitial' and a 'centered 
vacancy' equal Ge = 0.8 0.9 and Gtr = 0.85 0.95, (see Fig. 12(a)). It should be noted that in spite of prolonged 
anneahng of the system during 5 x 10^ MC steps at a temperature Ti = 0.007560To, which is just below melting, we 
did not find more complex isomers than point defects and quartets of dischnations. 

At the temperature T3 = 0.007565To, which is the melting temperature, (see Fig. 1(c)) we find pairs of dischnations 
which are created by the dissociation of dischnations which were bound into a quartet (see Fig. 5(b) of Ref. |39|). Pairs 
of dischnations can be found in an energy range of (0.216 ^ 0.285)fcsTo (group 2 in Fig. 12(a)). These defects destroy 
the translational order of the system, but conserve the bond-angular order. The energy range of (0.29 0.62)fcBTo 
(group 2) refers to multiple disclination pairs. Within this range of energy the bond-angular order factor remains 
larger than the translational one. For the energy range (0.62 l)kBTo (group 3) there are pairs of dischnations 
aggregating into chains, which are nothing else then grain boundaries (see Fig. 10) which have the effect of rotating 
one part of the crystal with respect to another. In this case, both order factors become very small. These defects occur 
in the system at intermediate temperature T = T2 when the system exhibits the remarkable property of switching 
between a quasi-liquid state and a crystalline state during our numerical simulation. These defects change the energy 
substantially, but allows the system to come back to a quasi-ordered state (see Fig. 11). At temperature T > O.OO8T0 
both order factors become small and the system transits to an isotropic fluid. 



VIII. TOPOLOGY OF DEFECTS IN THE SQUARE BILAYER CRYSTAL. 

In order to understand why the square lattice bilayer crystal has a considerable larger melting temperature as 
compared to e.g. the hexagonal lattice we investigated the various isomers which are created in the temperature 
interval Ti < T < T2 (see Fig. 2(c) for the location of the temperatures). For bilayer crystals the crystal structure 
and the topology of the defects is viewed as being composed by the top and the bottom staggered layers. Here, we 
took as an example v = 0.4 which has the largest melting temperature, i.e. T„iei = 0.01078ro- Note, that the energy 
of the defects which occurs in the square lattice depends on the interlayer distance. The scenario of defect formation 
is the same as in the case of a single layer crystal. First, for T < Ti = 0.01076 point defects are formed. These are 
the 'twisted bond' which is shown in Fig. 14(a) and has the energy E — 0.23fcBTo, a 'twisted triangle' and a 'twisted 
square' which are shown in Fig. 14(b,c), and have an energy E = 0.24fcBTo and E = 0.26kBTo, respectively (group 1 
in Fig. 12(b)). Each of these defects can appear separately. 

The next energy group includes the 'vacancy' and the 'interstitial' point defects which appear in pairs, see Fig. 14(d,f ) 
(group 1). The energy of a pair composed of a 'vacancy' and an 'interstitial' is E — 0.316fcB2o. These defects are 
followed by the formation of pairs of extended dislocations (see Fig. 5(e) of Ref. |^^). As seen from Fig. 12(b) these 
defects do not change substantially cither the bond-angular or the translational order factors. At T = r2 = 0.01078To 
the picture of defect structure changes crucially. Uncorrelated extended dislocations with non-zero Burgers vector 
and unbounded disclination pairs (see Fig. 5(e) of Ref. [^) are formed which causes a substantial decrease of the 
translational order. These isomers refer to the next energy range up to i? ~ lAksTo (group 2). Then with increasing 
density of dislocations and the appearance of single dischnations the system looses order and transits to an isotropic 
fluid (group 3) in Fig. 12(b). 

Fig. 12 clearly illustrates that the energy of the different defects in both lattice structures is substantially different. 
The defects which are able to destroy the translational and orientational order in the square lattice have a larger 
energy than corresponding one in the single layer. Even point defects in the square lattice crystal have a higher 
energy. Thus, the square type bilayer crystal requires larger energies in order to create defects which are responsible 
for the loss of the bond-orientational and the translational order. 



IX. CONCLUSION 



Monte-Carlo simulations of the melting transition of single and bilayer crystals were performed. The solid-liquid 
phase diagram for the five lattice structures which are stable in the bilayer system was constructed. We found an 
unexpected larger melting temperature for the square lattice. A comparison of the topology of the defects occurring 
in a hexagonal single layer crystal with these in the square bilayer crystal clarifies the enhanced stability of the square 
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lattice. In the case of the square lattice, all defects have a larger energy and consequently larger thermal energy is 
required to create them. 

An analysis of the bond-angular and translational order factors of a 2D crystal for the Coulomb, screened Coulomb, 

Lcnnard- Jones and l/r^"^ repulsive inter-particle interaction potentials allowed us to propose a new universal criterion 
for the melting temperature: melting occurs when the bond-angular correlation factor attains the value 0$ = 0.45. 
Using this criterion we showed that the melting temperature, TmeU can be obtained with sufficient accuracy within 
the harmonic approximation. 

It is shown that the potential energy of the system changes substantially within a very narrow range of temperatures 
around the melting temperature. We found the equilibrium state of the system nearby the melting point where the 
system alternates between the crystalline and the liquid state. The system is able to switch from the liquid to the 
solid state and back during the simulation time. Such a behaviour indicates that we reached the equilibrium state of 
the melting transition and hysteresis is absent for the melting of the Coulomb crystal at least for our finite size (448 
particles) system with periodic boundary conditions. 
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FIGURES 

FIG. 1. The relative mean square displacement of the particles (a), the height of the first peak of the pair correlation 
function (b), the potential energy (c) and the bond-angular (solid circles) and the translational (open circles) order 
factors (d) as function of temperature for the interlayer distance = 0. 

FIG. 2. The same as Fig. 1, but now for the interlayer distance v = 0.4, where the system forms a square lattice. 

FIG. 3. The same as Fig. 1, but now for the interlayer distance u = 0.8, where the system is in a staggered triangle 
lattice. 

FIG. 4. The pair correlation function of the square lattice bilayer crystal (v = 0.4) at T = O.OOOSTq (solid curve), 
T = 0.0095To (dotted curve) just before melting and T = 0.01 ITq (dashed curve) after melting where the system is 
in the liquid state. 

FIG. 5. The phase diagram of the bilayer Coulomb crystal for different screening lengths A = (open squares), 
A = 1 (solid circles), and A = 3 (open triangles). The curves are guides to the eye. The crystal structures are shown 
in the inserts where open (solid) symbols arc for the particles in the top (bottom) layer. 

FIG. 6. The structural phase diagram of the bilayer screened Coulomb crystal at T = as function of the screening 
length A. 

FIG. 7. The potential energy {E{T) — E{T = 0))/kBTo (a) and the bond angular order factor (b) as function of 
temperature for the interlayer distances: v = (solid squares), i/ = 0.4 (open triangles), for the screened Coulomb 
inter-particle potential with screening length A = 1. 

FIG. 8. The translational (a) and the bond-angular order factors (b) as function of temperature for the interlayer 
distances z/ = (A = 1-solid squares; A = 3-open circles), and f = 0.4 (A = 1-solid triangles; A = 3-open triangles), 
for the Lennard- Jones potential (solid rhombics) and for the repulsive potential (open rhombics). 

FIG. 9. The melting temperature in the harmonic approximation (solid curves) and from MC simulations (symbols): 
(a) for the square bilayer crystal as function of the interlayer distance for screening lengths A = (circles), A = 1 
(squares), A = 3 (triangles), and (b) for the hexagonal single layer crystal as function of the screening length. 

FIG. 10. Chain of dislocations (grain boundary): (a) in the vertical direction, and (b) in the horizontal direction 
in a single layer hexagonal lattice. 

FIG. 11. The bond-angular (squares) and the translational (open circles) order factors (a) and {E{T) — E{T = 
0))/fcs7o (b) as function of 'MC time' steps. 
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FIG. 12. The bond-angular (squares) and the translational order (circles) factors of the different defects in (a) a 
single layer crystal 1^ = 0, and (b) in the square lattice bilayer system for i' = 0.4. 

FIG. 13. Defects in a single layer crystal: (a) 'centred vacancy' and (b) 'centred interstitial'. 

FIG. 14. The point defects: (a) 'twisted bond', (b) 'twisted triangular', (c) 'twisted square', (d) 'vacancy', and (e) 
'interstitial' in the square lattice bilayer crystal. 
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